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ABSTRACT: 

The linearized barotropic primitive equations are used to investigate 
the phase speed errors which arise from space truncation. Phase speed 
errors are compared for second and fourth order approximations to the 
space derivatives. These errors are also investigated on two staggered 
grids . 
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1 . Introduction 



The objective of this report is to examine methods for improving the 
phase speeds of meteorological waves which are predicted by primitive 
equation models. In particular, the linearized phase speeds for second 
order and fourth order space difference approximations will be compared. 
When variables are computed at different grid points, computer time 
savings can be achieved. The phase speeds for two staggered grid 
schemes will be examined in an effort to find the scheme which will give 
the best phase prediction for the least computer time. 

2 . Phase Error Analysis 

The analysis will be applied to the barotropic primitive equations; the 
results can then be easily extended to the baroclinic primitive equations. 
The equation of motion and continuity equation can be written: 

-|t V + V . W + v 0 + fk x v = 0, (1) 

T'f + V . ( 0 V ) = 0 , (2) 

where 0 = gH; the quantity H is the height of the free surface. 

These equations can be linearized by writing: 

V = T U + u (x,t) 1 i + v (x,t) j , (3) 

L J ~ ~ 

0 = $- fUy + (p(x,t) , 

where U and 4- are constants and u, v and 0 are very small. If we 
substitute (3) into (1) and (2) and drop products of small quantities, we 
obtain the following set of linearized equations: 
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(6) 



We have neglected a term in the continuity equation which is proportional 
to 3 0/3y. The elimination of this term allows a simple separation of 
the modes, but it does not change their essential character. 

If we eliminate between Eqs . (4) - (6) we obtain the following equation 
for co : 
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The phase speeds of the waves can be obtained if we insert the wave 
ik. (x— ct) 

form , to = Ae into Eq . (7) . The three solutions are given by 
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<£ + f 2 /k 2 . 



(8) 



The first solution is the meteorological mode and the other two are 
inertial gravity waves . 

Our first objective is to compare the phase speeds of the meteorolog- 
ical modes computed from prediction equations which have second and 
fourth order approximations to the space derivatives. This analysis can 
be conveniently carried out in a semi-discrete system where the spatial 
variation is discrete and the time variation is continuous. The phase 
speeds computed from the semi-discrete system will be very close to 
those computed from the fully discrete system. It can be shown that for 
the meteorological mode, the difference between the phase speeds in the 
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discrete and the semi-discrete systems is proportional to (UAt/Ax) 3 . This 
ratio is very small because At must be small enough to keep the rapidly 
moving gravity waves stable. A stability analysis will be carried out 
later which will compare the values of A t which are needed for stability 
for the second and fourth order space difference schemes. 

If we write x = jAx and introduce second order differences into Eqs. 

(4) - ( 6 ), they become: 

iJij + u ( u j+i " u j-i ) + (<Pj+i ~ ) _ fv = 0, (9) 

a t i 

2 Ax 2 Ax 

5 V J + U + 1 ~ + fuj = 0, (10) 

9 t 2Ax 



1 <r j + u +i ~ cp j - 1 ) + 3 ( u i + 1 " u i- 1 ) = o . 
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3 t 2 A x 2 Ax 

Each field is now expressed in terms of a wave of wave number k: 

A , x ikAxj 
Uj = u (t) e 

A . . ikAxj 
Vj = v (t) e 

A ikAxj 
<Pj - <p (t) e 



(12) 



If we insert these expressions into (9) - (11) and divide out the exponential, 
we obtain: 
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where a = £inkAx)/Ax. Eliminate between these equations to obtain 



A 



the following equation for <p : 
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(16) 



(ti r - + iaU ) [ ("If +iaU ) +f a +a 2 $ ] £ = 



The solution is of the form 

A A -ikct 
u = u n e 



where 



aU/k 



c = 



(17) 



aU/k + V f* 



U/k + V f* / k 2 + ( a 2 / k 2 ) $ 

Since a -*k as A x -• 0, these expressions converge to the exact express- 
ions (8) . 

If we introduce fourth order space differencing into Eqs . (4) - (6), 



they become: 

Su < U 
+ 



( U J +1 - U 3-l) - ( U J + 2-U 3 -2 ) ] 



at i2Ax 

+ 12Ax [ ® ^ ) + i -( P i -i ) “ ( O j+s - <Pj-2 ) J _ l v j = 0 » 

[8 (Vj+i-V^!) - (Vj+2-Vj _2 ) ]+ fUj = 0, 
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+ f Uj = 0, 


(19) 
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{( ® J +1 <0 j -1 ) ” (<P J +-2 O ] -2 ) 1 

8 (Uj+i-u^ ) - (u j+3 -u 5 _ 2 ) 1 = 0. 



( 20 ) 



12 Ax 

If we introduce the relations (12) into these equations, they will take the 
same form as equations (13) - (15) if we replace a by y- 4/3 (sin kAx/Ax) 



- — [sin (2kAx)/Ax ] . In this case, the solutions for the phase speed 



4 



become: 



'(y/k) U 



C = 



(yk))U ± J f 2 /k 2 +(y 2 /k 2 ) * (21) 

These solutions also converge to the true phase speeds as Ax- 0. 

We are concerned with the phase speed of the meteorological mode 
whose exact value is U in this system. If U'is the phase speed of this 
mode in the semi-discrete system, we can write: 

( 22 ) 



U = U sin kA x/ (kAx) 

(second order space difference), 

l/ = U sin kAx/ (kAx) (1 + 2/3 sin 2 (kAx/2) ) 
(fourth order space difference). 



(23) 



Fig. 1 shows the errors incurred for each scheme as a function of kAx. 
This figure shows a considerable improvement for the fourth order scheme 
over the second order scheme except for the very short waves. 

Let us now compare the linear computational stability criteria for the 
second and fourth order schemes. Since the criteria is mainly determined 
by the gravity waves, we will only consider the following simplified 



system: 
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(24) 



( 25 ) 



Here we have set f = U = 0. If we write t = nAt and x = jAx and use 
second order differences, these equations become: 
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— n^n 



( ^J+ln - ^ j-ln ) , 



(26) 



u. j .1 u ' i 
jn+1 = jn-1 - 



<P 



jn+l = ^n-1 - TT *( Vln- U )-ln ) 



(27) 



We will employ the von Neumann method of stability analysis which is 

discussed by Haltiner (1971, p. 109). Introduce these Fourier relations: 
TT ikAxj -ikAxj 

u, = U e , (p. = P e . (2 8) 

Jn n J n n 

If we substitute these expressions into (26) and (27), we will obtain: 

Hv+l= U n-1- i sinkAx P , (29) 

Ax n 

n+1 = n- 1 - ^ ~ i sinkAxU . (30) 

Ax n 

Define the auxiliary variables : 



A = U , , B = P , 
n n-1 n n-I 



(31) 



and 



s = (2At/Ax) sinkAx, r = 2 $ (At/Ax) sinkAx . (32) 

With these definitions, the difference equations can be written in two- 
level form as the following matrix equation: 




The 4x4 matrix in this equation is called the amplification matrix. The 
eigenvalues of the amplification matrix satisfy the following equation: 
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X 4 - 2 ( 1 - sr/2 ) x 2 + 1 = 0 . 



(34) 



This equation can be solved for X 



V ( i - 



X = 1 - sr/2 ±V ( 1 - sr/2) -1 . (35) 

2 

When 4 i sr ^ 0, we find that I X 1=1. If sr > 4, it is clear that 
, 2 , 

IX I > 1 for one of the roots . The difference equations will satisfy the 
von Neumann necessary condition for stability if we can choose At and 



A x so that 4 a sr ^ 0 for all k. From (31) the product sr is 



sr 



= 4 4- ( ^ sin 2 kAx 



(36) 



It is clear from this expression that the second order scheme will be 
stable if 

At £ Ax/ $ ^ . (37) 

The difference forms of (24) and (25) with second order time differences 
and fourth order space differences are: 

U jn+1 U jn-1 6 Ax C^^j+lri ^j-ln^ ^j-i-2n ^j-2n^] , (38) 

^jn+l = * J »-1 ' £ ^ t 8( Vln- U j-ln * (U ) + 2n- U )-2n 1 ] . (39) 
If the Fourier relations (28) are substituted into these equations, we 



obtain: 

U n + 1 = U n-1- 3^ (8SinkAX ' Sln2kAx) p n ' (40) 

P n + 1 = P n-1 ‘ (8sinkAx - sin2kAx > U n . (41) 

The matrix equation (33) is the same for this case if we define 

s = At/(3Ax) ( 8 sinkAx - sin2kAx ) , 

\ (42) 

r = <fAt/(3Ax) (8 sinkAx - sin2kAx ) . 
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The solutions will be stable if 4 i sr i 0 for all k. If we introduce the 

relations (42), this criteria becomes: 

/ A t \ 2 / 4 1 \ 2 

4 ^ <1 J 4 ^ — sinkAx - — sin2kAx ) * 0 . (43) 

It can be shown that the maximum value of the quantity in the parenthesis 

is given by 

/ 4 1 \ 

( — sinkAx - — sin2kAx ) = 1.37 . (44) 

V 3 6 / max ' ' 

The stability criteria for the fourth order space difference scheme now 

becomes 

i 

At £ 0.73 Ax/4 8 . (45) 

If we compare (45) and (37), we see that the fourth order scheme requires 
a time step which is smaller than the time step required for the second 
order scheme. 

Let us now consider a mixture of second and fourth order differences 
which will allow a larger time step than the value given by (45) . If the 
derivation of the phase speeds (17) and (21) is followed through, it can be 
seen that the phase speed of the meteorological wave is independent of 
the difference approximation which is used for the pressure force term. 

The pressure force differencing affects only the gravity waves in this 
model. Consider a difference set in which the advection terms and the 
divergence term in the continuity equation are computed with fourth order 
differences and the pressure force term is computed with second order 
differences. With this set the phase speed of the meteorological wave 
will still be given by (21). We now apply this difference set to the 
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simplified Eqs . (24) and (25) in order to determine the critical value of the 
time step. The difference equations are given by (2 6 ) and (39) and the 
Fourier equations are given by (29) and (41). These equations can be 
written in the matrix form (33) if we define 

s = (2 At/ Ax ) sink Ax , (46) 

r = (<J>/3 ) (At/ Ax ) ( 8 sinkAx - sin2kAx ) . 

In this case, the stability criteria is 

4 s 4$(At/Ax ) 2 £ sinkAx (4/3 sinkAx - 1/6 sin2kAx) j ^0 , (47) 

for all k. It can be shown that the quantity in the brackets is non-negative 
and its maximum value is given by 

sinkAx ( 4/3 sinkAx - 7 - sin2kAx ) 1 = 1.35 . (48) 

6 A max 

The stability criteria now becomes 

i 

At ^0.86 Ax / <$ 2 . (49) 

This mixture of second and fourth order time differences allows a longer 
time step while giving fourth order accuracy in the phase speed of the 
meteorological wave . 

3 . Staggered Grid Schemes 

We will now examine the phase speed errors in two staggered grid 
schemes. The first, which is shown in Fig. 2, is used in the Arakawa- 
Mintz General Circulation Model (see Langlois and Kwok, 1969). 
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Fig. 2 

If we compare numerical predictions on this grid to predictions on a 
similar grid where all variables are carried at every grid point, we find a 
reduction in computer time and memory size of 4 to 1 . Although some 
averaging is required, the pressure force can be computed accurately 
since the differences are taken over a distance of 2 Ax. The advection 
terms are poorly represented, however, since the differences there are 
taken over a distance of 4 Ax. This gives a large error in the phase speed 
of the meteorological wave. The phase error for the meteorological wave 
of wavelength L on this grid is equal to the phase error of a wave of 
wavelength L/2 on the grid where all variables are carried at each point. 
Thus, for a wavelength of 6 Ax, the phase error ratio U'/U = 0.41 with 
this staggered grid while the value for the unstaggered scheme is U'/U = 
0.83. If fourth order differences are used with the staggered grid, the 
phase error ratio is improved to U'/U = 0.62. This is still poorer than 



10 



the value for the unstaggered grid. If we take into account the reduction 
in time step required by (45) or (49) and the increased computer time re- 
quired to compute the more complicated advection terms , the staggered 
scheme with fourth order difference will require less computer time than 
the unstaggered scheme. This staggered grid scheme is very useful for 
general circulation experiments or for ocean current prediction (see 
Haney (1971), but it is less useful for operational short range prediction 
where the minimization of phase errors is very important. 

Let us now consider the staggered grid array which is shown in Fig. 3 
which was examined by Winninghoff (1968). 
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Fig. 3 Ax 

This arrangement of grid points becomes the same as the arrangement 
shown in Fig. 2 if theaxesare rotated 45° and the grid size is increased. 
If we compare numerical predictions on this grid to predictions on a 
similar grid which is unstaggered, we find a reduction in computer time 
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and memory size of 2 to 1 . If a meteorological wave is propagating along 
either the x- or the y- axis on this grid, the phase error will be the same 
as on the unstaggered grid since the increments are taken over a distance 
of 2 Ax. This degree of accuracy requires the specific difference scheme 
which is presented later in this section. However, if the wave is propa- 
gating at an angle of 45° to one of the axes, the phase error will be 
larger since the differences are taken over a distance of 2 VTax. On the 
unstaggered grid, the errors are smaller since the differences are taken 
over a distance of VT Ax. Thus, for the 45° angle, the phase error on 
this staggered grid for a wave of wavelength L is equal to the phase error 
of a wave with wavelength y on the unstaggered grid. 

If we now introduce fourth order differences in the continuity equation 
and in the momentum advection terms, we will find an improvement in the 

phase speed of waves propagating along the axes as indicated in Fig. 1. 

o 

The propagation in the 45 directions will be improved by a similar 
amount, but it will still contain more error than with the unstaggered 
scheme. Even with fourth order space differences, this scheme will re- 
quire less computer time than for the unstaggered grid scheme. The 

phase propagation along the axes will be significantly better than for the 

o 

second order unstaggered scheme , while the propagation in the 45 
directions will be a little worse. Thus, this staggered scheme with fourth 
order space differences should give forecasts which are on the average 
better and require less computer time than the second order unstaggered 
scheme. Also, as pointed out by Winninghoff (1968), this staggered 
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scheme gives a better description of the geostrophic adjustment process 



than the unstaggered scheme. 

We will now present a set of difference equations for this staggered 
set of grid points . Equations (1) and (2) can be written in component 
form and combined to give the following three equations: 

J; (u0) + |- (uu0) + |- (vud) + <Z> - fv$ = 0 , (50) 

Ol O X 0 y O X 

+ lk ^ + "^y ^ + * + 0 ' ( 51 ) 

a <t> / 3 t + (u0) + ( v<t> ) = 0 . (52) 

o x a y 

For simplicity, the time variation will be left continuous. If we are to 
take full advantage of the placement of points on this grid, it is necessary 



to compute the flux divergences on a coordinate system which has been 
rotated 45° . Define the new coordinates as follows: 



II 


(x + y) , 


(53) 


n = -1 

/ 2 


(x - y) , 




and the corresponding velocity components are 




i = 1 

/2 


(u + v) , 


(54) 


b = 1 


(-U+v) . 





/2 



The flux divergence of the quantity s may be written 



-I- (u,«s) + £- (V4S) = iiiuj + ujmi , 

ox ay s £ or) 



( 55 ) 



where s = u or v in equations (50) and (51) . 
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Fig. 4 contains the new coordinates superimposed upon the staggered 
grid . 

Ax 

i i 

u,v 0 U,V 




j-1 ,m-l 



j,m-l 



j+1 ,m-l 



Fig. 4 



The flux terms (55) are computed in terms of fluxes across the sides of 
the diamond shown in Fig . 4 . In the finite difference formulation we 
require values for ^at some u, v points. Thus, we define 

6 . = 0 ■ , i + 0 . , + 0. , , + 0 . , j ( 

jm 4 \ j+lm j-lm jm+1 jm-1 1 

The finite difference forms of (50), (51), and (52) with second order 
differences are: 



14 



Qt + 8AX ^ 1 " + 1 +U J ”] [ (U * } 1+1* + ! + ( U ^ ) J » + ( V * ) J+ l. + l 

-[uj-i.-i +U JB ] [(u5) Jm + (U0 ) 2 . la - 1 +(v0) Jm + (viJj.j.-t ] 

+ [ U J - lm +1 +U JC ] [- ( U ^)j-lm + l - ( U ^) 3n + (V0)j_ ln+1 + (V0) JB ] 
Uj. + U J+lB _! ] [- (U0) JB - (U^ J+lffi _ x + M ) Jm + (v 0 ) j + in — i ] ^ 






+ 0 Jn ( 0 j +lm “ > / 2AX - f ( V 5 )j. = ( 5? ) 

5t" (V ^J" + V J"] [ (U0_) J+l^l + + ( V *) J+ lm + l + ( V ^) jm ] 

“ [ V J-i.-l + V JB ] [(U$) JB +(u5)j_ lB ^ + (V0) JB + (V0) j _ ln _ 1 ] 

+ [ V J-lm + l + V Jm ] [-iA0) 3 -l n+ l -(U0) JB + (V0)j- 1B+1 + (V0) JB ] 

- [ V jm + V 3 + l»-l ] [- ~ (ui) J+U -l + ^)j» + ( V ^) J+ l m -l ]^} 



+ 0 J» (Vi - ^--l) / 2Ax + f (U0) JB = 0, (58) 

iT" + Ttx [ (u ^3+i» " ( u 5) J - lm + (V0") JB+1 - (V0),.-! ] =0. (59) 

Fourth order differences can be introduced into any term if the differences 
are replaced by the weighted average of the differences computed on this 
grid and the differences computed on a grid with a double Ax. The weights 
are 4/3 for the small grid and -1/3 for the large grid. 

4 . Conclusions 

In this report we have shown that the use of fourth order differences in 
the momentum advection terms produces significant improvement in the 
phase speeds of short waves (see Fig. 1). The use of fourth order space 
differences, however, requires a time step which is 0.73 times the usual 
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time step to maintain computational stability. If the pressure force terms 
are computed with the second order differences and the other terms with 
fourth order differences, the time step can be increased to 0.86 times 
usual time step. This change in the differencing for the pressure force 
term does not affect the phase speed of the meteorological wave. 

Two staggered grid point arrangements were examined with respect to 
phase errors and computer time savings. The first scheme examined 
carries each variable at only one-quarter of the grid points with a corres- 
ponding savings in computer time. However, the advective phase errors 
for short waves are large in this system because the variables are differ- 
enced over twice as large a distance as with the unstaggered grid. Even 
the use of fourth order differencing will not bring the phase errors down to 
the value obtained from the unstaggered grid. The second scheme carries 
each variable at one-half of the grid points, which gives a savings of 
two in computer time. The phase errors on this grid are the same as on 
the unstaggered grid for waves moving along the coordinate axes. How- 
ever, for waves of wavelength L moving at an angle of 45° to the axes, 
the phase error in this staggered grid is equal to the error with a wave of 
wavelength yin the unstaggered grid. If fourth order differences are used, 
the overall phase propagation should be better than for the unstaggered 
grid. It is suggested that this scheme be tested with fourth order differ- 
ences in all terms except the pressure gradient term. This term should be 
computed with second order differences to allow a longer time step. The 
finite difference equations are given for a scheme which computes flux 
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divergence along lines which are at an angle of 45° to the coordinate 
axes. Winninghoff (1971) has tested a scheme on the sphere which has 
this grid point arrangement, but different advective approximations. 
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